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Abstract: It is currently believed that trabecular architecture 
in bone is result of adaptation process in which mechanical 
loads play a dominant role. Huiskes et al (2000) proposed a 
regulatory mechanism for this process. We implemented this 
model in a MATLAB code and studied effect of changing 
different parameters and other capabilities of such a 
mechanism. Starting from different initial configurations and 
by changing loading magnitude and direction in a plate model, 
we observed that this model always produces trabecular-like 
structures similar to those produced in simulation series 
performed by Huiskes and co-workers. In addition to 
confirming results of earlier simulations, such a general 
computer code can provide us with a better understanding of 
mechanisms and parameters involved in the adaptation 
process. Because 'stress shielding' is a major issue in long-term 
integrity of total hip replacements in the next simulation series 
we created a simple model of stress shielding around a 
prosthesis and investigated effect of changing bone cells related 
parameters. We also performed a comparison of this model and 
a topology optimization scheme. Surprisingly both procedures 
produced very similar results. Apparently this bone adaptation 
model has qualities similar to those of global optimization. 

INTRODUCTION 

Bone cells come into existence with the genetic blueprint 
for skeleton and sculpt it during growth until the skeletal 
design meets the loading demands. This process, termed 
bone adaptation, requires bone cells to detect mechanical 
signals and integrate these signals into proper changes in the 
bone architecture. Specialized cells, osteoclasts and 
osteoblasts, respectively are in charge of bone resorption and 
formation. 

Although mechanisms involved in the regulation of these 
'actor' cells are still vague, it is evident that mechanical 
feedback must be involved ([1]; [2]). 

By closely following the latest developments in bone 
physiology, many researchers have tried to develop 
mathematical models for the bone adaptation process. In 
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1976, Cowin and Hegedus [3] developed a sophisticated 
continuum, thermo -mechanical theory so-called adaptive 
elasticity theory. In this model, bone is defined as a porous 
medium with two phases: an elastic structure and an 
extracellular fluid. According to this model, adaptation is 
controlled by strain. Following adaptive elasticity theory, 
many other theories have been developed by others. For 
example, Rouhi et al. [4] replaced volume fraction by free 
surface density in the constitutive equations; Also Rouhi et 
al. [5] added a microcrack factor in the remodeling 
equations; Huiskes et al. [6] suggested that instead of strain, 
strain energy density (SED) can be used as a suitable 
mechanical stimulus for remodeling. 

In the beginning of the 21 st century, Huiskes and co- 
workers [7] developed a more refined semi-mechanistic 
theory. The new theory was based on the regulation scheme 
depicted in Figure 1. The purpose of this research is to 
simulate this model by developing a MATLAB code as the 
framework for further investigation of its capabilities and to 
have a better insight into adaptation of bone predicted by 
this model. In addition, effects of the parameter values in the 
model, Finite Element mesh and the external load and 
loading direction on the predicted morphology will be 
studied. Also two other simulation series including stress 
shielding simulation and a comparison with a topology 
optimization code will be performed. 
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Fig. 1 . The regulatory process proposed by Huiskes and co-workers [7] . 



METHODS 

In computational model presented by Huiskes et al. [7] 
bone tissue is assumed to contain » osteocyte cells per cubic 
millimeters located in the mineralized matrix with a total of 
A T in the domain considered. 

Each osteocyte /' measures a mechanical signal R(t) (in 
Jm^s 4 ), the strain energy density (rate), in its location. In 
turn the osteocyte recruits osteoblasts to form new bone 
depending on the difference between the measured signal 
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and a reference signal, k. The influence of an osteocyte on its 
environment is assumed to decrease exponentially with 
increasing distance from the osteoblasts. The influence of 
osteocytes /' on the osteoblast at location x is described by 
the spatial influence function 

f i (x) = exp(-d i (x)/D), (1) 
where d.(x) is the distance between osteocyte /' and location 
x. The parameter D represents the distance from an osteocyte 
at which location its effect has reduced to e~ l ' ; i.e. 36.8% 
[8]. 

The osteoblast recruitment stimulus is given by the 
stimulus value P(x,t) to which all osteocytes contribute 
relative to their distance from x, hence 

P(x,t) = f i f l (x) Mi R l (t) P) 

i 

where „ is mechanosensitivity of osteocyte i. 

The regulation of the relative density m in location x is 
governed by 
dm 
dt 
dm 
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= z{P(x,t) -k}- r oc for P (x ,t) > k 
■ for P(x,t)<,k (3) 



dt 

where T is a constant regulating the rate of the process, k is 

the threshold level for bone formation and r gc is the relative 

amount of mineral resorbed by osteoclasts per day [7] . 

This model includes a probability p of osteoclast 
activation per surface site at any time. This probability is 
assumed to be regulated either by the presence of 
microcracks or by disuse. The probability of resorption by 
microcracks was considered spatially random and was 
expressed as 

p(x,t) = constant, (4) 
where this constant was selected to be 10%. When assuming 
osteoclastic activation by disuse, the probability of 
resorption is higher in areas of lower strain. This strain 
dependent probability was formulated as 

p(x,t) = c[a-P(x,tj\ if P<a 

p(x,t) = if P>a (5) 
where c=12.5 and a=\.6 [7]. 

The elastic modulus E(x,f) at each location is calculated 
from density according to ([8]) 

E(x,t) = bxm(x,t) r (6) 
where b and y are constants. 

Simulations were performed on a 2x2 mm , two 
dimensional bone FE model. The structure was loaded by 
distributed stresses of 2 MP a at the edges and frequency of 
the load was sinusoidal at 1Hz [7]. 

We created the two-dimensional model by implementing 
above mathematical expressions in combination with FE 
formulation in a MATLAB code. 

The plate was meshed with 40x40 and 80x80 (only for the 
first simulation) four -node bilinear quadrilateral elements 
(Fig. 2). The sensor distribution was uniform associated with 
a sensor influence parameter of D = 0.01 mm and sensor 
density was assumed to be 1600/mm 2 [8]. In 
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Fig.2. Initial configuration (post-mineralized fetal stage) for (A) 40x40 
and (B) 80x80 element models. 

all simulations we only considered case of spatial random 
resorption with a probability of 10% because as stated by 
Huiskes et al [7] the only effect of strain-regulated 
resorption, is in the kinetics of the process. The homeostatic 
architecture is maintained but develops much faster in the 
simulation [7]. It means these two conditions are equivalent 
in the theory. 

The computations were performed on a Pentium 3 GHz 
CPU with 1 MB of RAM. 

RESULTS 

Ten different simulations were performed. The first and 
second ones tested whether the theory produces trabecular- 
like 2D configurations from conceptual initial architectures, 
representing bone in the post-mineralized fetal stage and 
from a uniform density (Fig. 3). Simulations were prolonged 
until no more gross architectural changes occurred, 
representing the homeostatic mature stage. 

In these two simulations the structures remodeled toward 
similar homeostatic configurations, in which trabecular-like 
structures were created and trabeculae were aligned to the 
loading direction (Fig. 3). Increasing the number of elements 
produced the same results. 

In the third simulation when the external load applied to 
the homeostatic architecture was rotated by 25°, the 
orientation of the trabeculae gradually reorientated as well, 
to align again with the external load (Fig. 4). 

In the forth simulation the regulatory mechanism was able 
to adapt the structure to alternative loading conditions. A 
25% increase in loads produced increased trabecular 
thickness and gain of bone mass (Fig. 5). 

The effect of overloading and unloading on trabecular 
adaptation was investigated in the fifth simulation where we 
artificially disconnected two trabeculae while the same 
externally applied load was maintained. The disconnected 
and therefore unloaded tarbeculae disappeared, while the 
neighboring overloaded trabeculae thickened (Fig. 6). 
In the sixth simulation by increasing Young's modulus for 
the upper side of square plate a simple model of stress- 
shielding of bone around prosthesis was simulated. Bone 
resorption happened close to the implant surface and this 
continued until a new homeostasis was obtained (Fig. 7). 

In the seventh, eighth and ninth simulations we changed 
bone resorption and formation related parameters. 
Decreasing r and p(x,t) by 50% had similar effects on 

bone loss due to stress shielding (Fig. 8 A, B). In both cases 
the severity of bone loss was reduced and trabeculae close 
to the implant were thickened. But in case of increasing 
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Fig.3. Transformation of morphology for (A) uniform density (B) post- 
mineralized fetal stage initial configurations 
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Fig.4. Orientation of the applied stresses was changed from i5°to 0° and 
the architecture adapted to align with the new stress orientation. 
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Fig.5. Effect of increasing loading magnitude by 25%. All trabeculae 
thicken. 
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FIG. 6. Two struts are artificially disconnected, while the external stress 
is maintained. After adaptation, the existing architecture is again adapted to 
the applied stresses by removal of the unloaded trabeculae and thickening of 
the overloaded trabeculae. 
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FIG. 7. Stress shielding is simulated by increasing the young modulus of 
upper side of the square. 

bone formation rate by 50% the change in rate of bone loss 
was little (Fig. 8 C). A trabecula in top right corner of the 
first two models was maintained but in the bone 
formation enhanced model this trabecula disappeared. 
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FIG.8. Effect of stress shielding coupled with (A) Decreasing r by 
50% (B) Decreasing p(x,t) by 50% and (C) increasing bone formation 
rate by 50%. 

In the tenth simulation the question was how will this 
remodeling theory compare to a topology optimization 
model? We obtained very similar final configurations (Fig. 
9) by using our code (bone adaptation theory prediction) and 
the topology optimization code by ole Sigmund [9]. 

Although this simulation was performed for a simple 
loading condition (compressive ramp load) but the similarity 
between the two final configurations is interesting. 
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Fig.9. (A) Two-dimensional plate model of bone tissue subjected to a 
compressive load as indicated (B) Architecture after applying the bone 
adaptation theory (C) Architecture after applying the topology optimization 
code. 

DISCUSSION AND CONCLUSION 

From first five simulations it was concluded that the 
theory was able to explain both modeling as well as 
remodeling of trabecular-like architectures as governed by 
external forces. 

Although we implemented the mathematical model in a 
computer code from scratch, we obtained final 
configurations very similar to the ones obtained by Huiskes 
et al. The differences between the figures shown here and 
the figures from Huiskes et al. [7] are due to a difference in 
the time increments and time constant used in the 
formulation. Osteoblasts and osteoclasts were often assumed 
to work together in bone. In a recent paper Pogoda et al [10] 
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discuss that bone resorption is independent of bone 
formation. Studying transgenic mice, they induced a near- 
complete and reversible osteoblast ablation. In these 
animals, osteoblast ablation led to a complete arrest of bone 
formation accompanied by bone loss, thus illustrating that, 
in mice, the bone resorption function is independent of bone 
formation [10]. They also provide clinical evidence that the 
sympathetic regulation of bone does exist in humans and 
plays a clinically important role in diseases. 

A merit of this theory is that it is not in contrast with these 
latest findings in bone physiology because in this theory 
there is still a place for biochemical factors and central 
control of bone remodeling. This model does not imply that 
biochemical pathways or central control are irrelevant, but it 
does show that mechanical feedback can be a potent 
coupling factor for the relevant biochemical processes to 
take place [7]. 

Simulation of stress-shielding was another capability of 
the theory. A major problem threatening the long-term 
integrity of total hip replacement is the loss of proximal 
bone often found around non-cemented stems in the long 
term [6]. It is generally accepted that 'stress shielding' is the 
cause for this problem: after implantation of the prosthesis 
the surrounding bone is partially 'shielded' from load 
carrying and starts to resorb. 

Our simulations also suggest (very qualitatively) that 
anti-resoptive strategy can be more effective than increasing 
bone formation in treatment of bone loss around implants. 

Another interesting result was that when the global 
optimization procedure was applied to structures similar to 
those we used in this research, they produced results that are 
analogous. The implication is that the trabecular architecture 
predicted by this bone adaptation model has qualities of 
global mechanical optimality. 

Maybe more complete theories will appear as our 
knowledge of cellular processes in body progresses which 
includes the biochemical pathways and central control of 
bone remodeling. But this only depends on the parallel 
advancing of many sciences including genetics, 
biomechanics, biochemistry and computational techniques. 
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